perm filename FAKE10.SAI[DOC,TOB] blob
sn#076270 filedate 1973-12-05 generic text, type T, neo UTF8
00100 BEGIN "EDGE_TEST"
00200 comment
00300 find best choice of conf and diff for Hueckel edge
00400 operator based on chosen false positive rate for
00500 no signal. For a given false positive rate, find
00600 choice which requires the least height step to
00700 achieve a 50% false negative rate. This uses
00800 steps generated by a Monte Carlo procedure;
00900
01000 DEFINE SAFE_ARRAY=1; comment this makes all arrays safe;
01100 REQUIRE "VISION.2[SYS,TOB]" SOURCE_FILE; comment vision language package;
01200 REQUIRE EDGE_FILE;
01300 REQUIRE EDGE_PRINT;
01400 REQUIRE TV_VARIABLES;
01500 REQUIRE "STEP4.SAI[FAK,TOB]" SOURCE_FILE;
01600
00100 comment
00200 conf=sqrt(s0sq/(s0sq+n0sq)) ie cosine
00300 diff=s0sq;
00400
00500 BOOLEAN SIMPLE PROCEDURE YES(STRING S);
00600 BEGIN "YES"
00700 S←ASK_STRING(S);
00800 IF EQU(S,"Y")∨EQU(S,"y")∨EQU(S,"YES")∨EQU(S,"yes")
00900 THEN RETURN (1) ELSE RETURN (0);
01000 END "YES";
01100 DEFINE NBHD=8; comment the size of the array with fake intensity values;
01200 LABEL FIRST_LOOP; ! entry for false positives and others;
01300 LABEL SECOND_LOOP; ! entry for false negatives;
01400 LABEL ERR_TEST; ! entry for angular and position errors;
01500 LABEL FALSE_NEG; ! entry for false negatives;
01600 LABEL AVG_DIFF; ! entry for avg value of diff;
01700 LABEL HISTO_LOOP; ! loop entry for histograms repeat with change scale;
01800 DEFINE MAX_COUNT=2000; ! the max number of test samples generated;
01900 INTEGER_ARRAY INTENSITY[1:NBHD*NBHD]; ! store fake intensities here;
02000 REAL_ARRAY N0SQA[1:MAX_COUNT+20]; ! store n0sq from edge oper here;
02100 REAL_ARRAY S0SQA[1:MAX_COUNT+20] ; ! store s0sq from edge oper here;
02200 INTEGER_ARRAY COUNT[0:100]; ! histogram of scaled s0sq values;
02300 SHORT INTEGER K,K1;
02400 INTEGER ENOUGH; ! how many samples (< max_count);
02500 SHORT REAL STEP_MEAN,STEP_AMPLITUDE,NOISE_STD;
02600 comment the mean,amplitude,noise of fake steps;
02700 SHORT REAL NSQ; ! scaling factor which depends on noise squared;
02800 SHORT REAL NOISE_SCALE; ! scaling factor for histograms;
02900 SHORT REAL THETA,RHO; ! angle,perp dist of edge from origin;
03000 SHORT INTEGER BYTES; ! number of bytes per word;
03100 SHORT REAL LAMBDA; ! wavelength (used to get phase near an edge);
03200 SHORT INTEGER N1,N5,KN1,KN5; comment 1% and 5% level counts and indices;
03300
03400 ENOUGH←MAX_COUNT;
03500
03600 comment now init the edge operator;
03700 CONF←.05;DIFF←.05; ! choose low values to reject nothing;
03800 BITS←6;
03900 MODIFY(BITS);
04000 LINLEN←WORDS(NBHD,BITS);
04100 TVWORD←GIOWD(INTENSITY);
04200 LSIDE←1;RSIDE←NBHD;
04300 FLINE←1;LLINE←NBHD;
04400 EJINIT;
04500 DIFF←0;CONF←.05;
04600
04700 comment set up values for fake edges with zero height;
04800 NOISE_STD←5.0; ! amplitude of gaussian noise;
04900 STEP_MEAN←(2↑BITS-1)*.5; ! offset to get positive intensities;
05000 STEP_AMPLITUDE←0; ! step height for fake edges;
05100 NOISE_SCALE←20.0;
05200
00100 IF EQU(ASK_STRING(" SKIP FALSE POSITIVES "),"Y") THEN GO TO SECOND_LOOP;
00200 FIRST_LOOP:
00300 MODIFY(NOISE_STD);
00400 MODIFY(ENOUGH); ! ask user if he wants another value;
00500 IF ENOUGH>MAX_COUNT THEN ENOUGH←MAX_COUNT;
00600 MODIFY(STEP_AMPLITUDE);
00700
00800 OUTSTR(CRLF&"FALSE POSITIVES ");
00900 FOR K1←1 STEP 20 UNTIL ENOUGH DO BEGIN "ACCUM1"
01000 FOR K←K1 STEP 1 UNTIL K1+20 DO BEGIN "ACCUM" ! report every 20 events;
01100 SHORT INTEGER N;
01200 LABEL CALLFK;
01300 CALLFK: FAKE(NBHD,NBHD,THETA,RHO,LAMBDA,
01400 STEP_MEAN,STEP_AMPLITUDE,NOISE_STD,NOISE,INTENSITY);
01500 comment generate an array with zero edge,random;
01600 EJLI(NBHD/2,NBHD/2,0,0); ! try the edge operator on it;
01700 N0SQA[K]←N0SQ; ! save the noise squared;
01800 S0SQA[K]←S0SQ; ! and save the signal squared;
01900 END "ACCUM";
02000 PRINT(K);
02100 END "ACCUM1";
00100 comment for a range of conf values, find the corresponding 1% and 5%
00200 false positives (above threshold with no edge signal) for fake
00300 edges whose noise sq and signal sq have been saved;
00400 N1←ENOUGH%100;
00500 MODIFY(N1);
00600 TCRLF;PRINT(N1);
00700 HISTO_LOOP:
00800 MODIFY(NOISE_SCALE);
00900 NSQ←NOISE_SCALE/(NOISE_STD↑2);
01000 FOR CONF←.95 STEP -.05 UNTIL .40 DO BEGIN "CONF_LOOP" ! loop over conf;
01100 SHORT REAL CONFSQ;
01200 SHORT INTEGER N;
01300 CONFSQ←CONF↑2;
01400 comment now form histograms over scaled s0sq values;
01500 FOR K←0 STEP 1 UNTIL 100 DO COUNT[K]←0; ! zero the histogram;
01600 FOR K←1 STEP 1 UNTIL ENOUGH DO BEGIN "ACCUM"
01700 SHORT REAL S;
01800 S←S0SQA[K];
01900 IF S/(S+N0SQA[K])>CONFSQ THEN BEGIN ! if above conf;
02000 N←S*NSQ+.5;
02100 IF N>100 THEN N←100;
02200 COUNT[N]←COUNT[N]+1;
02300 END;
02400 END "ACCUM";
02500 OUTSTR(CRLF&" OUTPUT HISTOGRAM");
02600 PRINT(CONF);
02700 INCHWL;
02800 IFC 0 THENC
02900 comment print histogram;
03000 FOR K←0 STEP 1 UNTIL 100 DO BEGIN
03100 TCRLF;PRINT(K);OUTSTR(" "&CVS(COUNT[K]));
03200 END;
03300 ENDC;
03400
03500 comment now find the 1% and 5% levels of the histogram tails;
03600 N←0;
03700 KN1←100;
03800 FOR K←100 STEP -1 UNTIL 0 DO BEGIN
03900 N←N+COUNT[K];
04000 IF N<N1 THEN KN1←K;
04100 END;
04200 OUTSTR(CRLF&"DIFF CUTOFF LIMITS "&CVG(N1/ENOUGH));
04300 OUTSTR(CRLF&CVG(KN1/NSQ));
04400 OUTSTR(CRLF&"DIFF CUTOFF LIMITS IN UNITS OF SQUARE STD DEV");
04500 OUTSTR(CRLF&CVG(KN1/NOISE_SCALE));
04600 PRINT(KN1);
04700 END "CONF_LOOP";
04800 IF EQU(ASK_STRING(CRLF&"REPEAT HISTOGRAMS "),"Y") THEN GO TO HISTO_LOOP;
04900 IF EQU(ASK_STRING(CRLF&"REPEAT VALUES"),"Y") THEN GO TO FIRST_LOOP;
00100 SECOND_LOOP:
00200 IF EQU(ASK_STRING(" SKIP FALSE NEGATIVES "),"Y") THEN GO TO ERR_TEST;
00300 BEGIN "MIN_STEP"
00400 LABEL TEST_LOOP;
00500 LABEL MIN_REPEAT;
00600 REAL STEP_FRACTION;
00700 INTEGER KQ;
00800 INTEGER PRINT_FLAG;
00900 PRELOAD_WITH 1.250,1.350,1.350,1.400,1.400,1.400,1.400,1.400,1.400,1.400,
01000 1.400,1.400,1.400;
01100 OWN REAL ARRAY DIFF_CUTOFF[1:20]; ! cutoffs for diff;
01200 PRINT_FLAG←1;
01300 MODIFY(PRINT_FLAG);
01400 MODIFY(NOISE_STD);
01500 NSQ←NOISE_STD↑2;
01600 STEP_FRACTION←1.0;
01700 MODIFY(ENOUGH); ! ask user if he wants another value;
01705 MODIFY(THETA);
01710 MODIFY(RHO);
01800 IF ENOUGH>MAX_COUNT THEN ENOUGH←MAX_COUNT;
01900 MIN_REPEAT:
02000 MODIFY(STEP_FRACTION);
02100 STEP_AMPLITUDE←NOISE_STD*STEP_FRACTION;
02200 PRINT(STEP_AMPLITUDE);
02300
02400 OUTSTR(CRLF&"FALSE NEGATIVES ");
02500 FOR K1←1 STEP 20 UNTIL ENOUGH DO BEGIN "ACCUM1"
02600 FOR K←K1 STEP 1 UNTIL K1+20 DO BEGIN "ACCUM" ! report every 20 events;
02700 SHORT INTEGER N;
02800 LABEL CALLFK;
02900 CALLFK: FAKE(NBHD,NBHD,THETA,RHO,LAMBDA,
03000 STEP_MEAN,STEP_AMPLITUDE,NOISE_STD,NOISE,INTENSITY);
03100 comment generate an array with zero edge,random;
03200 EJLI(NBHD/2,NBHD/2,0,0); ! try the edge operator on it;
03300 N0SQA[K]←N0SQ; ! save the noise squared;
03400 S0SQA[K]←S0SQ; ! and save the signal squared;
03500 END "ACCUM";
03600 PRINT(K);
03700 END "ACCUM1";
00100 comment for a range of conf values, and corresponding diff values, find
00200 number of rejections of fake edges;
00300
00400 TEST_LOOP:
00500 KQ←0;
00600 FOR CONF←.95 STEP -.05 UNTIL .40 DO BEGIN "CONF_LOOP" ! loop over conf;
00700 SHORT REAL CONFSQ;
00800 SHORT INTEGER N;
00900 KQ←KQ+1;
01000 N1←0;
01100 CONFSQ←CONF↑2;
01200 FOR K←1 STEP 1 UNTIL ENOUGH DO BEGIN "ACCUM"
01300 SHORT REAL S;
01400 S←S0SQA[K];
01500 IF S/(S+N0SQA[K])>CONFSQ THEN BEGIN ! if above conf;
01600 IF S>DIFF_CUTOFF[KQ]*NSQ THEN N1←N1+1
01700 ELSE IF PRINT_FLAG THEN OUTSTR(CRLF&CVG(S));
01800 END;
01900 END "ACCUM";
02000 TCRLF;
02100 PRINT(CONF);OUTSTR(" "&CVG(DIFF_CUTOFF[KQ]*NSQ));
02200 OUTSTR(" "&CVG(DIFF_CUTOFF[KQ]));
02300 PRINT(ENOUGH)PRINT(N1);
02400 INCHWL;
02500 END "CONF_LOOP";
02600
02700 IF EQU(ASK_STRING(CRLF&"REPEAT HISTOGRAMS "),"Y") THEN GO TO TEST_LOOP;
02800 END "MIN_STEP";
02900 IF EQU(ASK_STRING(CRLF&"REPEAT VALUES"),"Y") THEN GO TO SECOND_LOOP;
00100 comment find the angular and position errors at threshold;
00200 ERR_TEST:
00300 IF YES("ANGULAR AND POSITION ERRORS ") THEN BEGIN "ERR_BLOCK"
00400 INTEGER K;
00500 INTEGER N;
00600 SHORT REAL ALPHA,SUM_ALPHA,SUMSQ_ALPHA;
00700 SHORT REAL DELTA,SUM_DELTA,SUMSQ_DELTA;
00800 MODIFY(STEP_AMPLITUDE);
00900 MODIFY(NOISE_STD);
01000 MODIFY(RHO);
01100 MODIFY(THETA);
01200 FAKE(NBHD,NBHD,THETA,RHO,LAMBDA,
01300 STEP_MEAN,STEP_AMPLITUDE,NOISE_STD,NOISE,INTENSITY);
01400 comment generate an array with zero edge,random;
01500 OP_APPLY(NBHD/2,NBHD/2); ! try the edge operator on it;
01600 SUM_ALPHA←SUMSQ_ALPHA←SUM_DELTA←SUMSQ_DELTA←0;
01700 N←0;
01800 FOR K←1 STEP 1 UNTIL ENOUGH DO BEGIN "ERR_LOOP"
01900 FAKE(NBHD,NBHD,THETA,RHO,LAMBDA,
02000 STEP_MEAN,STEP_AMPLITUDE,NOISE_STD,NOISE,INTENSITY);
02100 comment generate an array with zero edge,random;
02200 IF EJLI(NBHD/2,NBHD/2,0,0)∧¬BCOMP
02300 THEN N←N+1; ! try the edge operator on it;
02400 IF K MOD 20 = 0 THEN PRINT(K);
02500 ALPHA←ATAN2(CX,CY);
02600 SUM_ALPHA←SUM_ALPHA+ALPHA;
02700 SUMSQ_ALPHA←SUMSQ_ALPHA+ALPHA*ALPHA;
02800 DELTA←SQRT(OPX↑2+OPY↑2)-RHO;
02900 SUM_DELTA←SUM_DELTA+DELTA;
03000 SUMSQ_DELTA←SUMSQ_DELTA+DELTA*DELTA;
03100 END "ERR_LOOP";
03200 STD_ALPHA←SQRT((SUMSQ_ALPHA-SUM_ALPHA↑2)/N);
03300 STD_DELTA←SQRT((SUMSQ_DELTA-SUM_DELTA↑2)/N);
03400 TCRLF;
03500 PRINT(STD_ALPHA);
03600 PRINT(STD_DELTA);
03700 END "ERR_BLOCK";
03800 IF YES("REPEAT ERROR BLOCK") THEN GO TO ERR_TEST;
00100 comment find the mean and std dev of s0sq/(noise_std↑2);
00200
00300 AVG_DIFF:
00400 BEGIN "AVG_SIGNAL"
00500 INTEGER K;
00600 SHORT REAL SUM_S0SQ,SQ_S0SQ;
00700 SHORT REAL STD_S0SQ; ! std dev of s0sq;
00800 OUTSTR(CRLF&"MEAN VALUE OF DIFF AND STD DEV");
00900 SUM_S0SQ←SQ_S0SQ←0;
01000 TCRLF;
01100 MODIFY(STEP_AMPLITUDE);
01200 MODIFY(ENOUGH);
01205 MODIFY(THETA);MODIFY(RHO);
01210 MODIFY(NOISE_STD);
01300 FOR K←1 STEP 20 UNTIL ENOUGH DO BEGIN "AVG_LOOP"
01400 INTEGER K1;
01500 PRINT(K);
01600 FOR K1←1 STEP 1 UNTIL 20 DO BEGIN "A1_LOOP"
01700 FAKE(NBHD,NBHD,THETA,RHO,LAMBDA,
01800 STEP_MEAN,STEP_AMPLITUDE,NOISE_STD,NOISE,INTENSITY);
01900 comment generate an array with edge;
02000 EJLI(NBHD/2,NBHD/2,0,0); ! try the edge operator on it;
02100 SUM_S0SQ←SUM_S0SQ+S0SQ;
02200 SQ_S0SQ←SQ_S0SQ+S0SQ*S0SQ;
02300 END "A1_LOOP";
02400 END "AVG_LOOP";
02500 STD_S0SQ←SQRT((SQ_S0SQ-SUM_S0SQ)/ENOUGH);
02600 OUTSTR(CRLF&"AVERAGE VALUE OF DIFF/SQUARED AMPLITUDE"&CRLF);
02700 OUTSTR(CVG(SUM_S0SQ/(ENOUGH*STEP_AMPLITUDE↑2)));
02800 OUTSTR(CRLF&"AVERAGE VALUE OF STD DEV OF DIFF/SQUARED AMPLITUDE "&CRLF);
02900 OUTSTR(CVG(STD_S0SQ/(STEP_AMPLITUDE↑2)));
03000 NOISE_STD←0;
03100 FAKE(NBHD,NBHD,THETA,RHO,LAMBDA,
03200 STEP_MEAN,STEP_AMPLITUDE,NOISE_STD,NOISE,INTENSITY);
03300 comment generate an array with edge;
03400 EJLI(NBHD/2,NBHD/2,0,0); ! try the edge operator on it;
03500 PRINT(NOISE_STD);PRINT(S0SQ);PRINT(STEP_AMPLITUDE);
03600 END "AVG_SIGNAL";
00100 comment now ask user which values he wants for diff and conf
00200 start with low values of step height, low enough that fewer
00300 than 50% of the fake edges are detected. If too few, then
00400 double step height. If too many, we have bracketed step
00500 height. We can thus halve the step and use a binary search
00600 to zero in on the right value;
00700
00800 FALSE_NEG:
00900 BEGIN "EDGE_SEARCH" comment search for edge height with 50% detection;
01000 REAL SEARCH_CONTROL; comment either 2 or 1/2,
01100 start below, double value,bracket,binary search;
01200 REAL STEP0,DSTEP;
01300 LABEL CHANGE,EXIT,CHANGE1;
01400 INTEGER K,J; ! loop count variables;
01500
01600
01700 comment ask for initial values;
01800 TCRLF;TCRLF;
01900 STEP0←0;MODIFY(STEP0);
02000 MODIFY(NOISE_STD);
02100 NOISE_SCALE←1.0;
02200 MODIFY(NOISE_SCALE);
02300 DIFF←(NOISE_STD↑2)*NOISE_SCALE;
02400 PRINT(DIFF);
02500 CONF←0.90;MODIFY(CONF);
02600
02700 STEP_AMPLITUDE←STEP0; COMMENT CHECK FALSE NEGATIVES;
02800 SEARCH_CONTROL←2.0;
02900 DSTEP←0.1*NOISE_STD;
03000 THETA←PI/2;
03100 RHO←0;
03200 CHANGE1:
03300 FOR J←1 STEP 1 UNTIL 10 DO BEGIN "TST1"
03400 INTEGER TOTAL, POSITIVE; ! count of total tries and positives;
03500 STEP_AMPLITUDE←STEP0+DSTEP;
03600 TCRLF;MODIFY(STEP_AMPLITUDE);
03700 MODIFY(DSTEP);
03800 POSITIVE←0;
03900 TOTAL←0;
04000 MODIFY(ENOUGH);
04100 FOR K←1 STEP 20 UNTIL ENOUGH DO BEGIN "TST2"
04200 INTEGER K1;
04300 PRINT(K);
04400 FOR K1←1 STEP 1 UNTIL 20 DO BEGIN "TST3"
04500 FAKE(NBHD,NBHD,THETA,RHO,LAMBDA,
04600 STEP_MEAN,STEP_AMPLITUDE,NOISE_STD,NOISY_STEP,INTENSITY);
04700 IF EJLI(NBHD/2,NBHD/2,0,0)∧¬BCOMP THEN POSITIVE←POSITIVE+1;
04800 TOTAL←TOTAL+1;
04900 END "TST3";
05000 END "TST2";
05100 PRINT(TOTAL);PRINT(POSITIVE);
05200 IF POSITIVE*2>TOTAL THEN SEARCH_CONTROL←0.5
05300 ELSE STEP0←STEP0+DSTEP;
05400 DSTEP←SEARCH_CONTROL*DSTEP;
05500 END "TST1";
05600
05700 CHANGE:
05800 MODIFY(THETA);
05900 MODIFY(RHO);
06000 LAMBDA←1.0;
06100 GO TO CHANGE1;
06200 OUTSTR("ANOTHER VALUE-TYPE Y"&CRLF);
06300 IF INCHWL≠"Y" THEN GO TO EXIT;
06400 MODIFY(RHO);
06500 MODIFY(THETA);
06600 GO TO CHANGE;
06700 EXIT:
06800 END "EDGE_SEARCH";
06900 END "EDGE_TEST";